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The propagation of a laser beam through turbulent media is modeled as a 
fractional Brownian motion (fBm). Time series corresponding to the center 
position of the laser spot (coordinates x and y) after traveling across air in 
turbulent motion, with different strength, are analyzed by the wavelet theory. 
Two quantifiers are calculated, the Hurst exponent, H, and the mean Nor- 
malized Total Wavelet Entropy, Swt ■ It is verified that both quantifiers give 
complementary information about the turbulence state. 

1. Introduction 

Wavelets-based tools have been shown to be well-suited to fractal processes and their 
analysis. This is mainly due to the fact that the wavelet transform incorporates in 
its definition two basic features, time and scale, which are of primary importance for 
fractal processes. Another remarks concerns the structure of the wavelet transform 
which, by construction, builds a signal by successive refinements, starting from a 
coarse approximation and adding finer and finer details at each step. Such a procedure 
is of course reminiscent of basic fractal constructions, thus suggesting we should make 
use of wavelets to synthesize a fractal pro cess 

Since the earlier leading work of MandelbrolP^ there is a lot of evidence that 
several facets of fully developed turbulent flows are fractals. Most of the applications of 
fractals to turbulence have been devoted to the study of subsets of the region occupied 
by turbulent flow where a given property is satisfied. For example, the turbulent /non- 
turbulent interface, some constant property surfaces (such as the iso-velocity and 
iso-concentration surfaces) and the structure of spatial distribution of dissipation 
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rates have been characterized as fractals or multifractals and the dimension has been 

measured! 6 * 7 * 8 * Then, it is natural to apply wavelet tools for analyzing atmospheric 
t ur bulence.™MTl 

The fractional Brownian motion (fBm) was introduced by Mandelbrot and Van 
NessP as an example of a process which contains an infinite domain of dependence 
with the intention of explaining the results reported by Hurst in 195 1.^ fBm as a 
model for turbulence is not newl 12 * 11 ^^ In particular, the fBm can be introduced as 
an alternative model for the turbulent index of refraction, and it can be shown that 
these processes allow to reconstruct most of the index properties.^ 

This work reports on the basic characteristics of the centroid position of a laser 
spot after the light has traveled through turbulent media. It is analyzed as a stochastic 
process, within a mathematically defined theoretical model: the fractional Brownian 
motion. 

That is, given a component of the centroid position / as a time t function, the 
geometrical nature of the graph (t,f(t)) is studied. Afterwards, the wavelet theory 
is used to characterize this centroid. Two quantifiers are obtained: the Hurst expo- 
nent, H, and the mean Normalized Total Wavelet Entropy, Swt- Their behaviors are 
compared; the analysis shows they describe different properties of the turbulence. 

2. Fractional Brownian motion 

The fractional Brownian motion of exponent H (Hurst exponent), with < H < 1, 
is a zero-mean Gaussian process Bh(x) with iGl such that: 

a. B H (0) = 0. 
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b. The difference Bh{x + Ax) — Bh{x) has finite dimensional normal distribution. 

This means that their increments are stationary Gaussian processes whose variance 
is 

E [(B H (x + Ax) - B H (x)) 2 } = a 2 \Ax\ 2H , (1) 

where a is a parameter — E[ ■ ] is, of course, the average. The nonstationary character 
of the fBm is evidenced by its covariance function given by 

E [B H (x)B H (y)} = a l{\x\ 2H + \y\ 2H - \x - y\ 2H ) . (2) 

The covariance of future increments with past ones is: 

Ph(Ax) = E[(B H (x) - B H {x - Ax)) ■ {B H {x + Ax) - B H (x))\ 

= a 2 (2 2H - l -l)\Ax\ 2H . (3) 

Note that pjj is independet of x and the parameter H determine the correlation of 
the increments. 

For H = 1/2 the correlation of past and future increments vanishes for all x, as is 
required for an independent random process. However, for H ^ 1/2 one has pn{Ax) ^ 
0. For H > 1/2 this quantity is positive and the process is called persistent.^ In this 
case, an increasing trend in the past implies on the average an increasing trend in 
the future and, conversely, a decreasing trend in the past implies on the average a 
continued decrease in the future. For if < 1/2 the process is called antipersistent. 
Now, an increasing trend in the past implies on the average a decreasing trend in the 
future, while a decreasing trend in the past makes on the average an increasing trend 
in the future — see Fig. [TJ 
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K. Helland and Van AttcP were the first to study the Hurst exponent in grid 
generated turbulence by applying a rescaled-range analysis — an usual measure of long- 
term persistence in geophysical time series — to turbulence velocity records measured 
within a wind tunel. They showed that there are some deviations from H = 1/2 — 
referred as the Hurst phenomenon. 

As a nonstationary process, the fBm does not have a spectrum defined in the 
usual sense; however, it is possible to define an empirical power spectrum of the form: 
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SbM) = JJfW+T- (4) 
This equation is not a valid power spectrum in the theory of stationary processes since 
it is a nonintegrable function, but it could be considered as a generalized spectrum. 
Through this interpretation of Eq. (J3| a self-similarity relation can be shown for the 
fBm. That is for Bh{x) with H and a parameters, one has that a H Bn(x/a — b) have 
the same finite dimensional distributions for all a > and b. The fractal dimension 
of sample functions of these processes is D = 2 — H. 

3. Time-Frequency Analysis 
3.1. Wavelet transform 

First introduced by Dennis Gaborp^ wavelet analysis is a method which relies on the 
introduction of an appropriate basis and a characterization of the signal by the dis- 
tribution of amplitude in this basis. If the basis is required to be a proper orthogonal 
basis, any arbitrary function can be uniquely decomposed and the decomposition can 
be inverted) 20 * 21 * 22 * 2 ^ Wavelet analysis is a suitable tool for detecting and character- 



izing specific phenomena in time and frequency planes. 

The wavelet is a simple and quickly vanishing (compactly supported) oscillating 
function. Unlike sine and cosine of Fourier analysis, which are precisely localized in 
frequency but extend infinitely in time, wavelets are relatively localized in both time 
and frequency. Furthemore, wavelets are band-limited; they are composed of not one 
but a relatively limited range of several frequencies. 

A wavelet family ip a fi is the set of elementary functions generated by dilations 
and translations of a unique admissible mother wavelet tp(t): 

^(t) = |a|- 1/2 ^(^), (5) 

where a, b G R, a ^ are the scale and translation parameters respectively, and t is the 
time. As a increases, the wavelet becomes narrower. Thus, one have a unique analytic 
pattern and its replications at different scales and with variable time localization. 

The continuous wavelet transform (CWT) of a signal S(t) G L 2 (M) (the space 
of real square summable functions) is defined as the correlation between the function 
S(t) with the family wavelet tp a ^ for each a and b: 

(W^S) (a, b) = \a\~ 1/2 j°° S(t)r (t^j dt = (S, Vv & >- (6) 

For special election of the mother wavelet function ip(t) and for the discrete set of 
parameters, dj = 2~ J and bj t k = 2 _J 7c, with j,k G Z (the set of integers) the family 

^j,k{t) = 2 j/2 ip(2 j t - k) j,keZ, (7) 

constitutes an orthonormal basis of the Hilbert space L 2 (M) consisting of finite-energy 
signals. 

6 



The correlated decimated discrete wavelet transform (DWT) provides a nonre- 
dundant representation of the signal, and the values (S, tp a ,b) constitute the coeffi- 
cients in a wavelet series. These wavelet coefficients provide relevant information in 
a simple way and a direct estimation of local energies at the different scales. More- 
over, the information can be organized in a hierarchical scheme of nested subspaces 
called multiresolution analysis in L 2 (M). In the present work, we employ orthogonal 
cubic spline functions as mother wavelets. Among several alternatives, cubic spline 
functions are symmetric and combine in a suitable proportion smoothness with nu- 
merical advantages. They have become a recommendable tool for representing natural 
sie 

na lcplEIl_fi 

gure 121 shows the cubic spline wavelet function used in this work. 
In what follows, the signal is assumed to be given by the sampled values S = 
{x(n),n = 1, • • • , M}, corresponding to an uniform time grid with sampling time T s . 
For simplicity, the sampling rate is taken as T s = 1. If the decomposition is carried 
out over all resolutions levels the wavelet expansion is: 



where N = log 2 M and the wavelet coefficients Cj(k) can be interpreted as the local 
residual errors between successive signal approximations at scales j and j + 1, and 
Tj(t) is the residual signal at scale j. It contains the information of the signal S(t) 
corresponding to the frequencies 2?~ 1 uj s < \uo\ < 2?uj s with uj s the sampling frequency. 

3.2. Relative Wavelet Energy 

Since the family {ipj t k{t)} is an orthonormal basis for L 2 (R), the concept of energy is 
linked with the usual notions derived from Fourier's theory. The wavelet coefficients 
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are given by Cj(k) = (Sjipj^) and the corresponding asociated energy will be its 
square. The energy at each resolution j = — 1, ■ ■ • , — N, will be 

^• = ^E^)> ( 9 ) 

•? k 

where Nj represents the number of wavelet coefficients at resolution j. The total 
energy can be obtained in the fashion 

£tot = J2 £ r ( 10 ) 

j<0 

Finally, we define the normalized pj values, which represent the Relative Wavelet 
Energy (RWE) by 

p j = £ j /£ tot , (11) 

for the resolution levels j = —1, —2, • ■ ■ , — N. The pj yield, at different scales, the 
probability distribution for the energy. Clearly, YljPj = 1 anc ^ ^ ne distribution {pj} 
can be considered as a time-scale density that constitutes a suitable tool for detecting 
and characterizing specific phenomena in both the time and the frequency planes. 

3.3. Wavelet entropy 

The Shannon entropy^ gives a useful criterion for analyzing and comparing probabil- 
ity distribution. It provides a measure of the information contained in any distribu- 
tion. We define the Normalized Total Wavelet Entropy (NTWS)- 27 1 28 1 29 1 as 

S WT = -^2p r ]jip j /£r u , (12) 

where 

S max = lnN. (13) 
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The NTWS appears to be a measure of the degree of order- disorder of the signal. 
It provides useful information about the underlying dynamical process associated with 
the signal.^ Indeed, a very ordered process can be represented by a periodic mono- 
frequency signal (signal with a narrow band spectrum). A wavelet representation 
of such a signal will be resolved at one unique wavelet resolution level, i. e., all 
RWE will be (almost) zero except at the wavelet resolution level which includes the 
representative signal frequency. For this special level the RWE will be (in our chosen 
energy units) almost equal to one. As a consequence, the NTWS will acquire a very 
small, vanishing value. A signal generated by a totally random process or chaotic one 
can be taken as representative of a very disordered behavior. This kind of signal will 
have a wavelet representation with significant contributions coming from all frequency 
bands. Moreover, one could expect that all contributions will be of the same order. 
Consequently, the RWE will be almost equal at all resolutions levels, and the NTWS 
will acquire its maximum possible value. 

Figure El presents two different relative wavelet energy (probability) distribution 
corresponding to five wavelet resolution levels (j = —5, . . . , —1). It is clear from the 
figure that distribution A presents broad band spectrum. In contrast, distribution B 
shows a clear dominance of the resolution level j = —2. According to the description 
given above, for the NTWS the following relations can be expected: NTWS (A) > 
NTWS (B). This is observed in the numerical values given at the figure. 
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3.4. Wavelet quantifiers time evolution 

In order to follow the temporal evolution of the quantifiers defined above, RWE and 
NTWS, the analyzed signal is divided into % non-overlapping temporal windows with 
length L (i = 1, • • • , Nt, with N? = M/L). Afterwards, appropriate signal- values for 
these quantifiers are assigned to the middle point of each time window. In the case 
of a diadic wavelet decomposition, the number of wavelet coefficients at resolution 
level j is two times smaller than at the previous, j + 1, one. The minimum length of 
the temporal window L will therefore include at least one wavelet coefficient at each 
level. 

The wavelet energy at resolution level j for the time window i is given by 



3 fc=(i-l)-L+l 

where Nj represents the number of wavelet coefficients at resolution level j corre- 
sponding to the time window i; while the total energy in this time window will be 




i-L 




(14) 




(15) 



i<o 



The time evolution of RWE and NTWS will be given by: 



P 3 £j /£••, 



(16) 




(17) 



In order to obtain a quantifier for the whole time period under analysis^ the 
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temporal average is evaluated. The temporal average of NTWS is given by 

(S W t) = S WTi ( 18 ) 

T 1=1 

and for the wavelet energy at resolution level j 

N T 

1 i=i 

then the total wavelet energy temporal average is defined as 

<S*) = I>i>- (20) 
i<o 

In consequence, a mean probability distribution {qj} representative for the whole time 
interval (the complete signal) can be defined as 

qj = (£,)/(£ tot ), (21) 

with Ylj Qj — 1 an d the corresponding mean NTWS as 

s WT = -J2qr^ qj /s ma *. (22) 

j<0 

4. Fractional Brownian motion and Wavelet Transform 

A relevant property of the wavelet based multiresolution analysis is the stationary 
character of the wavelet coeficient series corresponding to each level j of resolution.^ 
Another important property is that the reconstruction of the original time series from 
the stationary series of wavelet coefficients reproduces the original signal with small 
error. 

In relation with fractional Brownian motion it can be shown that: 

a. fBm is nonstationary but the wavelet coefficients are stationary at each scale^^ 
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b. fBm exhibits positive long-range correlation in the range 1/2 < H < 1 but 
wavelet coefficients have a correlation which can be arbitrarily reduced; 

c. the self-similarity of fBm is directly reflected in its wavelet coefficients, whose 
variance varies as a power law as a function of scale ^SE3 

log 2 {E [C](k) \ Bh ] } oc -(2H + (23) 

All the above properties evidence that wavelet analysis is naturally well-suited to 
fBm. Each of them provides in fact a key ingredient for a problem of major importance 
when analyzing fBm: the estimation of the Hurst exponent or of the related spectral 
exponent a = 2H + 1.EE3 Starting from the above mentioned observations that, in 
the wavelet transform, variance progression follows a power law across scales 

V[j]=E[q(k) \ BH ] ~2"* (24) 

one can simply make use of the empirical variance estimators at scale j (based on 
Nj = 2~- 7 M coefficients for a sample of total length M) 

Nj 

V\3] = ^Y J C 2 J {k)\ BH =£,\ BH , (25) 
^ fc=i 

This is made possible because the fBm, when decomposed via the wavelet transform, 
becomes stationary at each scale. 

In this way an estimator of the parameter H can be obtained by: a) estimating 
the variance of the wavelet coefficients with Eq. (|23j) : b) plotting log 2 {V r [j]} versus j 
and fitting a minimum square line. The slope of the line give the estimator of H. 

Wavelet-based estimators dedicated to fBm can be viewed as versatile general- 
izations of previous techniques. An important feature of fBm is that its increments 
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are stationary and such that the structure function is proportional to |Ax| — see 
Eq. (0) — thus suggesting that this variance should be estimated in order to find H . 
Although feasible, this approach is faced with a difficulty due to long-range depen- 
dence. Classical (empirical) variance estimators are obviously poor estimators in such 
a context and specific estimators have to be designed F"^ 4 * 

5. Experimental Setup and Data Adquisition 

Time series corresponding to the fluctuations in the position of a laser beam's spot 
(wandering) over a screen, after propagation through laboratory generated convective 
turbulence, were recorded with a position sensitive detector located as screen at the 
end of the path. Twenty records of 5000 spot beam coordinates measurements every 
five minutes were obtained. Temperature along the laser beam path, for each record, 
were also measured and stored. 

The experimental measures were performed in the laboratory by producing con- 
vective turbulence over a path of lenght L = 1.5m with two electric heaters in a 
row and covering the path. We use a lOmW continuous wave He-Ne laser (Melles 
Griot Model 05-LHP-991), wavelenght laser beam A = 632. 8nm (red), with point- 
ing stability, after 15 minutes of warmup, less than 30/xm and a beam divergence of 
1.24±5%mrad. 

Each electric heater was 50cm long and 800W of power. The height of the laser 
beam's propagation path was lm above the electrical heaters and a box of expanded 
polyurethane was used as a thermal protection for the equipment. Three thermome- 
ters were allocated on the top of the box, 10cm above the laser beam. They were used 
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to determine the temperature in three positions along the propagation path, see Fig. 
01 The geometrical experimental arrangement was similar to that used by Consortini 
et al.' 1 "' and, Consortini and O'DonnellP^ 

The position sensor has a relative accuracy of the order of 2.5/xm so that very 
small position fluctuation can be measured. It was interfaced to a computer which 
allowed to measure at a rate of about 800 samples per second (M = 5000 laser spot 
coordinates in approximately 7s). Thus, with these coordinates stored on a hard disk 
the Hurst exponent and the mean NTWS were computed off line. 

Three different intensity levels of convective turbulence were generated by chang- 
ing the amount of heat dissipated for each electrical heater: 

a. Normal turbulence, the electrical heaters were off; 

b. Soft turbulence, each electrical heaters dissipated half of its available power; 

c. Hard turbulence, each electrical heaters dissipated at the maximum available 
power. 

Figure 03 shows the temperature for the three turbulence levels over the twenty 
records. Note that there are temperature gradients over the twenty records for the 
three turbulence levels and that the temperature (Tl, T2 and T3, see Fig. HJ) is 
not homogeneously distributed. So, an uniform flux of warm air and uniformity of 
turbulence along the path was not present. 
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6. Results and discussion 

We transformed the twenty time records corresponding to the laser beam x and y- 
coodinates (M = 5000), under the three turbulence conditions, to the time-frequency 
domain by means of orthogonal discrete wavelet transform (ODWT) 2 " * 21 * 22 * 2 ^ obtain- 
ing in this way the corresponding wavelet coefficients Cj(k) series. In the present work, 
we consider eight wavelet resolution levels and cubic spline mother wavelet EZEEES3 
Note that the wavelet coefficients were non-overlapping for each scale. 

The wavelet coefficients were squared to obtain the associted wavelet energy and 
the estimator of wavelet coefficient variance V[j}- To evaluate the Hurst exponent 
the procedure described at the end of Sec. |U is used. Figure |H1 presents log 2 {V A [j]} 
versus the resolution level j, for the (x, y)-coordinates and for the three turbulence 
conditions. The vertical lines represent the scaling region used for the square fitting. 
The values of the coefficient H obtained were averaged over the twenty records at 
each turbulence condition for the (x, ^-coordinates respectively. 

On the other hand, each signal under analysis was divided in time windows of 
length L = 256 and for each one the probability wavelet energy distribution was 
evaluated. With these values the mean wavelet energy distribution and the mean 
NTWS, Swt, were obtained. This quantifier is taken as representative of the order- 
disorder of the whole signal. As before, the obtained values Swt were averaged over 
the twenty records at each turbulence condition for the (x, ^-coordinates respectively. 

A comparison between average Hurst exponent and mean NTWS is given in Fig. 
[7| It is known the entropy is a measure of the order of a given system. In this case, 
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the mean NTWS shows the same behavior for both coordinates. As the turbulence 
increases the mean NTWS does the same. But, the Hurst exponent discriminates 
between coordinates: for the x-axis no noticeable change is observed, and for the 
y-axis the Hurst exponent decreases (increasing the roughness, see Fig. with the 
increasing turbulence. 

Thus, it is observed that the Hurst exponent is sensitive to the mean flow of the 
warm air. It distinguishes the anisotropy characteristic of the convective turbulence. 
Because the entropy measures the order obviously it does not detect the mean flow. 

In any case, the value of the obtained quantities indicates that memory as well 
as self-similarity and scale invariance are significant property of these time series. 

In relation with the Hurst exponent a useful generalization consists in allowing the 
singularity exponent to become time-dependent H(t), thus generating a new process 
B H (t)(x) such that 

E [{B H{t) {t(x + Ax) - B m {t(x)) 2 ] = a 2 \Ax\ 2H{t) . (26) 

In such a situation the increment process is no longer stationary. Then, it is impos- 
sible to globally apply the technique mentioned above for estimating H(t). Provided 
that variations of H(t) are smooth enough, the time series can be divided in time 
windows where this requirement is satisfied and for each one the similar techniques 
based on time-scale energy distributions can be applied locally.^ In any case, we must 
emphasize that for each one an scaling region must be defined. The time evolution 
of NTWS could be easily implemented EBEa Moreover, the NTWS is capable of de- 
tecting changes in a nonstationary signal due to the localization characteristic of the 
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wavelet transform, the computational time is significantly shorter since the algorithm 
involves the use of wavelet transform in a multiresolution framework and, the NTWS 
is parameter-free and scaling region is not necessary for its evaluation. 

In the next future the intention of our project is to make experiments varying 
the intensity of turbulence by adjusting the voltage applied to the heating element 
using a voltage controller, in order to study the quantifiers temporal evolution and 
characterize in a quantitative way the dynamics of the process. In another hand, the 
same study will be doing in outdoor experiments at different moments of the day to 
characterize the ground atmospheric turbulence. 
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Fig. 1. Sample paths of fractional Brownian motions with H = 0.3 (antipersis- 
tent), H = 0.5 (standard brownian motion), and H = 0.7 (persistent). These 
graphs were obtained by using the software FRACLABp^ 
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Fig. 2. Cubic spline wavelet. 
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Fig. 3. Relative wavelet energy (probability) distributions correspond- 
ing to five wavelet resolution levels (j = — 5,...,— 1). Distribu- 
tion A, { Pj } = {0.05,0.10,0.30,0.35,0.20}; distribution B, { Pj } = 
{0.03,0.10,0.12,0.70,0.05}. The NTWS values for this distribution are 
NTWS(A) = 0.888 and NTWS(E) = 0.614. 
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Fig. 4. Experimental setup. 
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Fig. 5. Temperatures in positions 1, 2 and 3 for hard (a), soft (b) and normal 
(c) turbulence. 
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Fig. 6. Evaluation of Hurst exponent for laser beam x-coordinate (left column) 
and y-coordinate, for hard (top), soft (center) and normal (bottom) turbulence. 
The vertical lines represent the scaling region used in the evaluation of H. 
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Fig. 7. Averages for Hurst exponents and mean NTWS (with their respective 
standard deviation errors) for the x and y coordinates in the case of normal, 
soft, and hard turbulence. 
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